#include "ellipsincludes.h"
#include "fresnel.h"
#include "film.h"

// 29/5/2010 add menu driver to routines
//
typedef struct ri_d {
	double refractive_index;
	complex double thickness;
} ri_d_t;


// menu choices
char *method_choices[] = {
	"Fixed n, find thickness",
	"Fit refractive index and thickness",
	"View/Change defaults",
	"Previous menu"
};
WINDOW *create_newwin(int height, int width, int starty, int startx)
{	WINDOW *local_win;
	
	local_win = newwin(height, width, starty, startx);
	box(local_win, 0 , 0);		/* 0, 0 gives default characters 
	 * for the vertical and horizontal
	 * lines			*/
	wrefresh(local_win);		/* Show that box 		*/
	
	return local_win;
}


int CalcD2(double angle, double delta, double tanp, double complex n1, double complex n2, double complex n3, double lambda, double complex *d2a, double complex *d2b)
// calculate the thickness of layer 2 using McCrackin's method.

// angle in radians
// n1, n2, n3 refractive indices
// lambda wavelength

// need to make angles (and Snell() complex)
// use cos^2=1-sin^2?


{
//	double  a1,a2,a3; // angles in the media 
	double complex cos1,cos2,cos3; // cosines of angles in the media
	double complex rp12,rp23,rs12,rs23; // fresnel terms
	double complex A, B, C; // quadratic coefficients
	double complex X1,X2;

//	a1=cos(angle);
//	a2=snell(n1,n2,a1);
//	a3=snell(n2,n3,a2);
//	printw("psi=%g    delta=%g\n\n",atan(tanp)*180./M_PI,delta*180./M_PI);
	cos1=ccos(angle);
	cos2=cos_snell(n1,n2,angle);
	cos3=cos_snell(n1,n3,angle);
	
	// now calculate the fresnel terms
	// rp
	rp12=rp(n1,n2,cos1,cos2);
//	printf("rp12=%g +I%g\n",creal(rp12),cimag(rp12));
	rp23=rp(n2,n3,cos2,cos3);
	
	//rs
	rs12=rs(n1,n2,cos1,cos2);
	rs23=rs(n2,n3,cos2,cos3);
	
	// quadratic terms
	
//	A=-rp12*rp23*rs23*cexp(delta*I)*tanp + rp23*rs12*rs23;
	A=tanp*cexp(delta*I)*rp12*rp23*rs23-rp23*rs12*rs23;
	
//	B=-rp12*rp23*rs12*cexp(delta*I)*tanp + rp12*rs12*rs23 - rs23*cexp(delta*I)*tanp + rp23;
	B=tanp*cexp(delta*I)*(rs23+rp12*rp23*rs12)-(rp12*rs12*rs23+rp23);
//	C=-rs12*cexp(delta*I)*tanp + rp12;
	C=tanp*cexp(delta*I)*rs12 - rp12;
	
	
	X1=(-B+csqrt((B*B)-(4.*A*C)))/(2.*A);
#ifdef _debugging_
	printw("X1=%g %g\n",X1);
	printw("log X1=%g %g\n",clog(X1));
	printw("abs X1=%g\n",cabs(X1));
#endif		   
	X2=(-B-csqrt((B*B)-(4.*A*C)))/(2.*A);
//	printf("X1=%g +I*%g\n",creal(X1),cimag(X1));
	*d2a=((lambda*I)/(4.*M_PI*n2*cos2))*clog(X1);
#ifdef _debugging_
	printw("Re{d2a=%g\n",creal(*d2a));
#endif
	// d2a=((lambda*I*clog(X1))/(4.*M_PI*n2*cos(a2)));
	*d2b=((lambda*I)/(4.*M_PI*n2*cos2))*clog(X2);
#ifdef _debugging_
	printw("thick=%g\n",(lambda/(4.*M_PI*n2*cos2))*creal(I*clog(X1)));
#endif
//	printf("d2a=%g +I*%g\n",creal(((lambda*I*clog(X1))/(4.*M_PI*n2*cos(a2)))),cimag(((lambda*I*clog(X1))/(4.*M_PI*n2*cos(a2)))));
 return 0;
}

int getpsidelta_solve(WINDOW *mainwnd, ellipspar_t *ellipsparams)
{
	// prompts user for psi delta values then calls solvepsidelta to solve for n and thickness
	double psi;
	double delta;
	//double angle;
	char texta[20];// user input angle
	psi=0.0;
	delta=0.0;
	mvprintw(1,1,"Input film psi and delta (both in degrees)");
	mvprintw(3,1,"Psi?");
	
	refresh();
	echo();
	getstr(texta);
	psi=atof(texta);
	
	mvprintw(5,1,"Delta?:");
	
	refresh();
	echo();
	getstr(texta);
	delta=atof(texta);
	
	
	mvprintw(7,1,"Psi=%g, Delta=%g\n",psi,delta);
	refresh();
	//convert to radians
	
	delta=M_PI*delta/180.;
	psi=M_PI*psi/180.;

		
	
	solvepsidelta(mainwnd,ellipsparams,psi,delta);
	
	return 0;
}

int bisect(double angle, double delta,double  tanp,double n1,double nmin,double nmax, double  n3,double lambda,double complex * d2a, double complex * d2b, double* nresult)
{
	// returns zero if all went well uses pointers to get values back to calling function
	// returns 1 if there was a problem with the Calc2D routine
	// returns 100 if the is no solution in the range nmin to nmax
	// returns 200 if there is no convergence.
	int error;
	int finished =1;
	double complex best1, best2, best3;
	while (finished<30)
	{
		// bisection style loop to find best parameters
		if((error=CalcD2(angle,  delta,  tanp,  n1, nmin,  n3,  lambda,  d2a, d2b))!=0) 
		{
			printw("Calc2D problem (at nmin)");
			refresh();
			return 1;
		}
		// choose best result at nmin
		best1=*d2a;
		if (abs(cimag(*d2a))>abs(cimag(*d2b))) best1 = *d2b;
		
		if((error=CalcD2(angle,  delta,  tanp,  n1, nmax,  n3,  lambda,  d2a, d2b))!=0)
		{
			printw("Calc2D problem (at nmax)");
			refresh();
			return 1;
		}
		//choose best result at nmax
		best2=*d2a;
		if ( abs(cimag(*d2a)) > abs(cimag(*d2b)) ) best2 = *d2b;
		
		// test to see if there is a sign change 
		if (cimag(best1)*cimag(best2)>0.)
		{
			// no sign change therefore stop running
			//			printw("No solution found\n");
			//			printw("1. thickness %g  RI %g      2. thickness %g  RI %g\n",creal(best1), creal(nmin), creal(best2),creal(nmax));
			//			refresh();
			return 100;
		}
		// must have a sign change to have got this far
		
		
		
		
		//	printw("%d 1. thickness %g  RI %g      2. thickness %g  RI %g\n",finished,creal(best1), creal(nmin), creal(best2),creal(nmax));
		//	refresh();
		// look for convergence	
		if ((abs(cimag(best1))<0.05)||(abs(cimag(best2))<0.05))
		{
			//have converged
			// 
			printw("Result found\n");
			printw("%d. thickness %0.0f  RI %0.2f      2. thickness %0.0f  RI %0.2f\n",finished,creal(best1), creal(nmin), creal(best2),creal(nmax));
			printw("abs(cimag(best1)=%g  abs(cimag(best2)=%g\n",abs(cimag(best1)),abs(cimag(best2)));
			refresh();
			*nresult=(nmin+nmax)/2.;
			return(0);
		}
		// No convergence so work on new range
		// calculate median refractive index
		double test=(nmin + nmax)/2.0;
		
		if((error=CalcD2(angle,  delta,  tanp,  n1, test,  n3,  lambda,  d2a, d2b))!=0) return 1;
		
		
		
		best3=*d2a;
		if (abs(cimag(*d2a))>abs(cimag(*d2b))) best3 = *d2b;
		
		// choose range to search.
		if(cimag(best1)*cimag(best3)>0.)
		{
			// answer lies between median and nmax
			nmin = test;
		} else
		{
			nmax=test;
		}
		
		finished++;
	}
	return 200;	
}

int solvepsidelta(WINDOW *mainwnd,ellipspar_t *ellipsparams, double psi, double delta)
{
	// solves psi delta pair for n and thickness
	// psi and delta in radians
	//
	
	
	double complex d2a,d2b;
	double lambda; // wavelength in angstromms
	double tanp=(tan(psi));
	double complex n1,n2,n3; // refractive indices
	double n2r=1.46; // real part of film refractive index
	double complex best1;
	double nmin,nmax,nresult; // range of refractive indices
	double angle; // ellipsometer angle

	int error;
	scrollok(mainwnd, TRUE);
	lambda=ellipsparams->ellipspar_values[5];
	n3=ellipsparams->ellipspar_values[1] + ellipsparams->ellipspar_values[2]*I; // substrate values
	
	n1=ellipsparams->ellipspar_values[3] +ellipsparams->ellipspar_values[4]*I; // air
	
	n2=1.46 - 0.*I; //start values (SiO2)
	angle=M_PI*ellipsparams->ellipspar_values[0]/180.;
	nmin=0.;
	nmax=0.;
	nresult=0.;
	tanp=tan(psi);
while (1) {
	
	
	//find out how user wants to solve the system
	int method=MethodMenu(mainwnd);
	// have choice, now operate on it:
	switch (method)
	{
		case 1:
		// "Fixed n, find thickness"
			if((error=CalcD2(angle,  delta,  tanp,  n1, n2,  n3,  lambda,  &d2a, &d2b))!=0) return 1;
#ifdef _debugging_
			printw("d2a=%g + %gi\n", creal(d2a), cimag(d2a));
			printw("d2b=%g + %gi\n", creal(d2b), cimag(d2b));
#endif
			refresh();
			best1=d2a;
			if (abs(cimag(d2a))>abs(cimag(d2b))) best1 = d2b;
			
			printw("the best solution is: t=%g A, with an error factor of %g\n", creal(best1), cimag(best1));
			refresh();
			printw("Press any other key to return to menu.");
			refresh();
			
			int current_getch = getch();
			
			
			break;
		case 2:
			// vary n and d
			// try varying n to find fit
			// Try to bracket the solution
			get_ri_limits(mainwnd, &nmin, &nmax);
			printw("nmin=%g  nmax=%g\n",nmin,nmax);
			refresh();
			
			// scan through possible values of refractve index and find the one with the most physically reasonable thickness value	
			search(angle, delta, tanp, ellipsparams->ellipspar_values[3]+I*ellipsparams->ellipspar_values[4], nmin, nmax, ellipsparams->ellipspar_values[1]+I*ellipsparams->ellipspar_values[2], ellipsparams->ellipspar_values[5], &d2a, &nresult);
			wprintw(mainwnd, "Best value found: thickness %f, refractive index: %f, confidence parameter: %f",creal(d2a),nresult,cimag(d2a));
			wprintw(mainwnd,"press any key to return to menu");
			getch();
			
			break;
		case 3:
			//change fixed refractive index
			getDoubleValue(mainwnd,"Film refractive index",&n2r,MIN_RI,MAX_RI);
			n2=n2r+0.*I;
			break;
			case 4:
			return 0;
				break;

	}
}
	// set up values
#ifdef _debugging_
	printw("ri subs=%g + %gi\n", creal(n3), cimag(n3));
	printw("ri top=%g + %gi\n", creal(n1), cimag(n1));
	printw("wavelength %g\n",lambda);
	printw("Angle=%g    %f\n",angle,ellipsparams->ellipspar_values[0]);
#endif
// assume layer is SiO2



return (0);
}

int search(double angle, double delta,double  tanp,double complex n1,
		   double  nmin,double nmax, double complex n3,double lambda,
		   double complex * d2, double *nresult)
{
	
	double nstep=0.005;
	double best = HUGE_VAL;
	
	double complex d2a=0.;
	double complex d2b=0.;
	ri_d_t best_pair; //stores best refractive index / thickness pair
	
	best_pair.thickness=0.;
	best_pair.refractive_index=0.;
	
	for (double nval=nmin; nval<=nmax;nval+=nstep) 
	{
		
		CalcD2(angle, delta, tanp, n1, nval, n3, lambda, &d2a, &d2b);
		//find the best of d2a&b
		// see if this best solution is better than the current global best.
		if (cabs(cimag(d2a))<=best) 
			{
			// found a new best value
			best_pair.refractive_index=nval;
			best_pair.thickness=d2a;
			best=cabs(cimag(d2a));
			//	printw("nval:%f       %f\n",nval,cabs(cimag(d2a)));
			}
		if (cabs(cimag(d2b))<best)
					  {
						  best_pair.refractive_index=nval;
						  best_pair.thickness=d2b;
						  best=cabs(cimag(d2b));
						 // printw("nval:%f\n",nval);
					  }
		refresh();
	}
									
	// now have best value								
	// 
	*d2=best_pair.thickness;
	*nresult=best_pair.refractive_index;
									
									
	return (0);
					
}



int MethodMenu(WINDOW *mainwnd)
//int main()
{	int startx=0;
	int starty=0;
	int n_choices = sizeof(method_choices) / sizeof(char *);
	
	WINDOW *method_menu_win;
	int highlight = 1;
	int choice = 0;
	int c;
	
	initscr();
	clear();
	noecho();
	cbreak();	/* Line buffering disabled. pass on everything */
	startx = (80 - WIDTH) / 2;
	starty = (24 - HEIGHT) / 2;
	
	method_menu_win = newwin(HEIGHT, WIDTH, starty, startx);
	keypad(method_menu_win, TRUE);
	mvprintw(0, 0, "Use arrow keys to go up and down, Press enter to select a choice");
	refresh();
	print_choice_menu(method_menu_win, highlight,n_choices);
	while(1)
	{	c = wgetch(method_menu_win);
		switch(c)
		{	case 65:
				if(highlight == 1)
					highlight = n_choices;
				else
					--highlight;
				break;
			case 66:
				if(highlight == n_choices)
					highlight = 1;
				else 
					++highlight;
				break;
			case 10:
				choice = highlight;
				break;
			default:
				refresh();
				break;
		}
		print_choice_menu(method_menu_win, highlight,n_choices);
		if(choice != 0)	/* User chose to come out of the infinite loop */
		{
			clrtoeol();
			refresh();
			delwin(method_menu_win);
			return choice;
			break;
		}	
		
		
	}
}
	void print_choice_menu(WINDOW *method_menu_win, int highlight,int n_choices)
	{
		int x, y, i;	
		
		x = 2;
		y = 2;
		box(method_menu_win, 0, 0);
		for(i = 0; i < n_choices; ++i)
		{	if(highlight == i + 1) // Highlight the current choice 
		{	wattron(method_menu_win, A_REVERSE); 
			mvwprintw(method_menu_win, y, x, "%s", method_choices[i]);
			wattroff(method_menu_win, A_REVERSE);
		}
		else
			mvwprintw(method_menu_win, y, x, "%s", method_choices[i]);
			++y;
		}
		wrefresh(method_menu_win);
		return;
	}
